Spatial correlation reverses the compound effect of multiple stressors on rocky shore biofilm

Abstract Understanding how multifactorial fluctuating environments affect species and communities remains one of the major challenges in ecology. The spatial configuration of the environment is known to generate complex patterns of correlation among multiple stressors. However, to what extent the spatial correlation between simultaneously fluctuating variables affects ecological assemblages in real‐world conditions remains poorly understood. Here, we use field experiments and simulations to assess the influence of spatial correlation of two relevant climate variables – warming and sediment deposition following heavy precipitation – on the biomass and photosynthetic activity of rocky intertidal biofilm. First, we used a response‐surface design experiment to establish the relation between biofilm, warming, and sediment deposition in the field. Second, we used the response surface to generate predictions of biofilm performance under different scenarios of warming and sediment correlation. Finally, we tested the predicted outcomes by manipulating the degree of correlation between the two climate variables in a second field experiment. Simulations stemming from the experimentally derived response surface showed how the degree and direction (positive or negative) of spatial correlation between warming and sediment deposition ultimately determined the nonlinear response of biofilm biomass (but not photosynthetic activity) to fluctuating levels of the two climate variables. Experimental results corroborated these predictions, probing the buffering effect of negative spatial correlation against extreme levels of warming and sediment deposition. Together, these results indicate that consideration of nonlinear response functions and local‐scale patterns of correlation between climate drivers can improve our understanding and ability to predict ecological responses to multiple processes in heterogeneous environments.


| INTRODUC TI ON
Ecosystems face multiple interacting natural and human-driven disturbances whose occurrence, magnitude, and impact vary in space and time (Gunderson et al., 2016;Jentsch et al., 2009;Wernberg et al., 2013). Predicting the ecological impacts of environmental fluctuations is essential as anthropogenic climate change alters the frequency and intensity of climate extremes (Drijfhout et al., 2015;Kirtman et al., 2013). Most notably, climate variability is already impacting on gross primary production at continental scales, and it is altering patterns of biodiversity at local scales (Ciais et al., 2013;Dornelas et al., 2014;Franks et al., 2007). In the last two decades, laboratory and field experiments have documented strong ecological effects of environmental fluctuations on species and assemblages through changes in frequency, variance, and autocorrelation of disturbance, resource supply, or other variables (Benedetti-Cecchi, 2003;Benedetti-Cecchi et al., 2006;Bernhardt et al., 2018;Bertocci et al., 2005;Crain et al., 2008;Gunderson et al., 2016;Lawson et al., 2015). One general mechanism explaining these results is the prevalence of nonlinear response functions relating ecological and environmental variables. Nonlinear responses, where a change in the input can generate a disproportionate change in the output, are ubiquitous in nature (Denny, 2017;Zhang et al., 2015).
If f(x) is a nonlinear function of an environmental variable x, nonlinearity causes a mismatch between the expected response under average conditions, f(x), and the integrated response under variable condition, f(x), such that f(x) ≠ f(x). This mathematical property of nonlinear functions is known as Jensen's inequality or nonlinear averaging (Jensen 1906). The sign of the inequality is positive (negative) for accelerating (decelerating) response functions, whereas the magnitude of change depends on the degree of nonlinearity and the amount of variability in × (Chesson, 2012;Ruel & Ayres, 1999).
The great emphasis on large temporal and spatial scales of variation in environmental conditions has overshadowed the importance of local-scale environmental variability (Helmuth et al., 2006;Sears et al., 2011). Yet, predictions and ecological outcomes at local scales can be dramatically different from those generated by global climate models (Nadeau et al., 2017). Landscape configuration and heterogeneity have been shown to play an essential role in modulating the impact of large-scale climate forcing on both plants and animals (Dong et al., 2017;Lehtilä et al., 2020;Ohler et al., 2020;Sunday et al., 2014). Spatial heterogeneity in the thermal environment may create local refugia, allowing species to survive during unusually harsh conditions (e.g., during heatwaves). For example, local-scale microclimatic conditions in plant communities can determine longterm resilience in rear edge forests again heatwaves and drought pressures (Carnicer et al., 2021). In addition, spatial heterogeneity may also ensure the persistence of thermally sensitive species (Angilletta, 2009). However, the extent to which spatially variable thermal regimes determine the performance of organisms in natural environments remains largely understudied (Dowd et al., 2015).
Despite the widely recognized importance of environmental variability in influencing organisms' performance, few studies have explored the consequences of Jensen's inequality in a multifactorial context, incorporating the effect of correlation between variables (Koussoroplis & Wacker, 2016;Koussoroplis et al., 2019;Pincebourde et al., 2012). The simultaneous effect of multiple factors may produce nonlinearities that would remain undetected in a unifactorial scenario (Koussoroplis et al., 2017). These nonlinearities may interact with environmental variance and the degree and direction of the correlation among environmental factors, leading to a deviation between organism's performance under constant and variable conditions. Direct evidence of the role of correlation in modulating organisms' performance stems from recent laboratory and mesocosm experiments (Koussoroplis & Wacker, 2016;Koussoroplis et al., 2019;Pincebourde et al., 2012). For instance, Koussoroplis and Wacker (2016)  Rocky intertidal habitats have been extensively used as model systems to test cornerstone hypotheses in ecology and to unravel the influence of multiple interacting processes (Hawkins et al., 2020). Plants and animals living on rocky shores are exposed to a mosaic of environmental conditions where key variables such as temperature and wave action often covary at experimentally tractable spatial scales (Helmuth et al., 2006;Lima & Wethey, 2012;Hawkins et al., 2020). The small size and short life span of many organisms living on rocky shores further facilitate the analysis of multiple stressors and their correlation in space or time. For example, recent studies have used epilithic microphytobenthos (hereafter, biofilm) to show how the temporal clustering of extreme events of warming and sediment deposition can promote legacy effects and drive populations to collapse . These and other studies have documented strong negative effects of warming and sediment accretion following heavy rains on rocky intertidal organisms (Harley, 2003;Kordas et al., 2015;Vaselli et al., 2008). Substratum complexity, generated by emergent surfaces mixed with heterogenous areas with depressions, cracks, and crevices in the rock, can result in different patterns of spatial correlation between aerial temperature and sediment. For example, positive correlation may occur on emergent rocks on sunny days and calm sea, whereas cracks and crevices where sediment is more persistent during rough sea conditions may introduce negative correlation.
We used a combination of experimental and simulation approaches to evaluate how spatial correlation between two important climate variables, warming and sediment accretion following run-off, modulated the performance of rocky intertidal biofilm. As a first step, we used a response surface design involving 16 combinations of warming and sediment deposition to derive a response surface (hereafter RS) relating the biomass and photosynthetic activity of biofilm to the two variables in the field. Then, we used the empirically derived RS to simulate the performance of biofilm under different correlation scenarios of warming and sediment deposition.
We expected the performance of biofilm to differ between constant and variable conditions owing to nonlinearities in the RS. Finally, we tested the predictions originating from our simulations in a second field experiment in which we manipulated the spatial correlation (positive or negative) and intensity of warming and sediment deposition in a factorial experiment. Our results provide the first empirical evidence of how spatial correlation between interacting stressors and nonlinear effects drives small-scale ecological responses of primary producers in real-world conditions and indicate that local-scale patterns between climate variables may play a crucial role in predicting ecological responses to multiple processes in heterogeneous environments.

| Study site
The study was done along the coast of Calafuria (Livorno, 43° 30′

| Experiment 1: Derivation of the response surface
We used a response-surface experimental design involving 16 combinations of warming and sediment deposition to build a warming-sediment response surface (RS) (Figure 1a). The experiment was conducted between May and August 2017. In May 2017, 48 plots consisting of patches of rock 40 × 40 cm fully covered by biofilm were marked at their corners with rawl plugs inserted into the rock for future relocation. Three replicated plots were randomly allocated to each combination of four levels of warming crossed with four levels of sediment deposition. Temperature and sediment were manipulated following methodologies developed in previous studies (Dal . The warming factor included a control (ambient temperature) and three levels of elevated temperatures (+5°C, +10°C, and +15°C above ambient levels). Plots were heated with aluminium chambers equipped with stoves ( Figure S1), with warming levels chosen to reflect a wide range of temperatures experienced by the biofilm at the study site, from common to rare.
We characterized individual warming levels as the averaged return time of positive thermal anomalies using a 66-year time series of temperature measurements (Appendix S1, Figure S2 Experimental treatments were applied only once as a single pulse; due to the impossibility of treating all the 51 plots (48 used to derive the response surface and the three CA plots) on the same day, sets of 5-6 randomly chosen plots were treated in each of 6 days within a month.
Biofilm biomass and photosynthetic activity were evaluated after 7, 14, and 21 days since the start of the experiment. Biomass was determined by means of an image-based remote sensing technique that uses chlorophyll a concentration as a proxy (μg chl a cm −2 ). Chlorophyll a was estimated from a ratio of reflectance at near-infrared (NIR) and red bands (Ratio Vegetational Index -RVI) by means of a IR-sensible camera (ADC), following the method proposed by Murphy and colleagues (Murphy et al., 2009). NIR/red ratios are related to chlorophyll a by a linear relationship, calculated on the basis of laboratory chlorophyll a extraction from Calafuria sandstone cores (Dal Bello et al., 2015). Each photo was then handled with a java-routine in ImageJ software to haphazardly select six subplots and provide a mean value of biofilm biomass for each plot (Schneider et al., 2012).
The physiological status of the photosynthetic apparatus of biofilm was assessed through a portable underwater pulse-amplitudemodulated (PAM) fluorometer (Diving-PAM, Walz). Maximum photochemical efficiency after 5′ of dark adaptation (henceforth dark yield) and effective quantum yield of photosystem II in actinic light (henceforth light yield) were used as a proxy of photosynthetic efficiency and stress, respectively. Within each experimental plot, three and six measurements were haphazardly taken for light and dark yield, respectively. Sampling had an average duration of about 2 h and started around 2.5 h after sunrise and ended at midday. Dark yield measurements were taken after 5 min of dark adaptation, while light yield was measured under natural light condition.
We used a generalized additive model (GAM) to derive the RS from the experimental data.
We obtained a single RS for each response variable (Chl a, dark and light yield) by taking averages over the three sampling dates.
Data were modeled as a function of three smoothers of nominal levels of warming (W) and sediment deposition (S) and their interactions (W × S) (with identity link and Gaussian error distribution): where the Y ijk is the value of the response variable (Chl a, dark and light yield) in replicate k, sediment level j, and warming level i, ß 0 is the intercept, s 1 and s 2 are thin plate regression splines describing the individual effect of warming and sediment, te is tensor product smooth term modeling the interaction between warming and sediment, and is the Gaussian error term. Smoother terms were selected through Flow chart illustrating the main steps of the study. (a) To derive a response surface, we performed a full factorial experiment crossing four warming intensities (+0°C, +5°C, +10°C, and +15°C above ambient air temperature, corresponding to 24.8°C, 30.4°C, and 40.0°C absolute temperatures, respectively) with four levels of sediment deposition (0 cm, +0.5 cm, +1.0 cm, and + 1.5 cm thick layers of sediment deployed over the plots). (b) Hypothetical warming and sediment response surface. The response surface was used to generate predictions of biofilm performance under variable (B(Warm. , Sed. )) and constant (B Warm., Sed. ) conditions for different scenarios of spatial correlation between warming and sediment deposition. The inset shows the effect of nonlinear averaging, -i.e., the deviation of biofilm performance due to the nonlinearity characterizing the warming-sediment RS -for a specific combination of warming and sediment deposition. Warming and sediment deposition are expected to follow a normal bivariate distribution. (c) We performed a field experiment to test the predicted outcomes of the simulations by manipulating the intensity and spatial correlation of warming and sediment deposition along experimental transects. (d) Transects consisted of three contiguous quadrats of 40 × 40 cm. Spatial correlation between warming and sediment deposition was generated by varying the levels of the two variables along the quadrats in a transect. The panel shows a high intensity treatment with an average level of warming of +15.5°C above ambient temperature, an average layer of sediment of 1.5 cm, and a level of spatial correlation between the two variables of −1. a generalized cross-validation procedure. The number of knots was set to 4, which corresponded to the number of levels of predictors variables. Model assumptions were assessed visually using plots of residuals vs. fitted values, box plots of residuals vs. experimental conditions, and QQ plots of standardized residuals vs. normal quantiles (Faraway, 2016). GAM fitting was performed using the function gam of package mgcv in R 3.5.1 (Wood et al., 2016).
Potential artifacts due to shading effects during warming sessions were assessed through a t-test contrasting control and CA plots.

| Simulating from the response surface
The RS experiment identified a significant nonlinear response of biofilm biomass to warming and sediment deposition (in interaction), but not for photosynthetic activity (see RESULTS -Experiment 1: Derivation of the response surface). Thus, we used the RS model derived for biomass to simulate the response of biofilm to changes in warming and sediment deposition under constant and variable scenarios and for different patterns of correlation between the two stressors.
These simulations allowed us to explore the interactive effects between nonlinearities and environmental variance on biofilm performance and the modulating effect of spatial correlation. We , where T and S corresponded to the chosen prediction values, defined the strength of correlation between predictors, and T and S were the standard deviations of temperature and sediment thickness estimated from field measurements during the experiment ( T = 4.54°C and S = 0.41 cm) ( Figure S4). In particular, T was calculated from temperatures in control plots, while S was calculated from data of the thickness of sediment deposits taken after heavy rainfall events. Thus, the two standard deviations T and S quantified the natural levels of spatial variability of warming and sediment deposition at the study site. We used the mean over 1000 simulations to obtain the predicted value of biofilm biomass for each combination of temperature and sediment deposition from the RS model. This procedure was repeated for different values of correlation ( ) ranging from −1 to 1, with increments of 0.2.
In addition to recording changes in biofilm biomass in the various scenarios, we used the simulated values to compute the total variance effect (TVE), a quantity that summarizes the effect of variance and correlation between stressors (Koussoroplis & Wacker, 2016). The

| Experiment 2: Testing predictions
Simulations generated quantitative predictions on the nonlinear response of biofilm biomass to changes in variance and correlation between warming and sediment deposition. We tested some of these predictions in a second experiment in which we manipulated the intensity and spatial correlation between stressors. In August 2019, we marked 12 transects each consisting of three contiguous quadrats (40 × 40 cm) in areas originally covered by biofilm ( Figure S1d).
The experiment had a factorial design with two levels of correlation (+1 and −1) crossed with two levels of intensity (low: +5°C of warming and 0.5 cm of sediment deposition; high: +15°C of warming and 1.5 cm of sediment deposition) and three replicate transects in each treatment combination. Each transect constituted a single experimental unit in which spatial correlation was manipulated by exposing each of the three quadrats to a specific combination of warming and sediment deposition. For example, in the positive correlation treatment, quadrats exposed to high (low) temperatures also received high (low) levels of sediment deposition. In contrast, in the negative correlation treatment, quadrats exposed to low (high) warming also received high (low) levels of sediment deposition (Figure 1d). Initial values of temperature and sediment for the three quadrats in each transect with designated levels of correlation (+1 and −1) were generated by sampling a bivariate normal distribution (Figures 1d, S1c,d).
Due to the small sample size involved (the three quadrats in a transect), this procedure often resulted in levels of spatial correlation that were lower than the nominal level. In these instances, the values of the two variables were adjusted arbitrarily to obtain the desired level of correlation. As in the first experiment, three additional transects were used as controls for artifacts (CA) to assess the potential effect of shading on biofilm during the heating sessions.
To compute the TVE, each of the four combinations of intensity and correlation of stressors was matched with an independent set of three replicate transects (12 in total: three replicates × two levels of correlation × two levels of intensity) with zero correlation between warming and sediment deposition. These additional transects provided the constant condition at the denominator of Equation 2 and (2) were obtained by imparting the same levels of warming and sediment deposition across the three quadrats in each transect. We recognize that although these treatments had zero nominal covariance, realized correlations could differ from zero owing to small variation in treatments levels across quadrats in a transect. Nevertheless, since Equation 2 uses the mean values of biofilm biomass across treatments, it is the realized mean correlation that must be close to zero.
We verified this expectation by recording the amount of warming and sediment deposition imparted to all quadrats during the experiment and computing the corresponding correlation values ( Figure S4).
Chl a was used as a surrogate for biofilm biomass and was mea-
The RS showed four distinct regions with contrasting nonlinear responses of biofilm biomass to warming and sediment deposition. A concave-8(0-10°C) and at intermediate sediment deposition (1 cm) and extreme warming (15°C) (Figure 2a). In contrast, a concave-up relation was evident at low to moderate levels of warming (0-10°C), at moderate to extreme levels of sediment deposition (1-1.5 cm) and, to a less extent, under extreme warming and low sediment deposition (Figure 2a). A response surface fitted to absolute temperatures ( Figure S5) provided a similar outcome, although it explained slightly less variation (delta temperatures in Figure 2a: Neither dark nor light yield exhibited a significant relationship with warming and sediment deposition ( Figure S6, Table S2). Shading effects or other artifacts due to the heating chambers were not detected for any of the three variables examined ( Figure S7, Table S4).

| Simulations
The strength and direction of correlation between warming and sediment deposition modulated the compounded effect of these stressors on biofilm biomass (Figure 3)

| Experiment 2: Testing predictions
The impact of warming and sediment deposition on biofilm was strongly modulated by the degree of spatial correlation between these stressors (Figure 5a). The analysis identified a significant interaction between correlation and treatment intensity (LMEM for the Correlation × Intensity interaction: t = 2.98, p < .05, df = 12; for the Correlation × Intensity interaction: t = 3.23, p < .01, df = 12; Figure 5b, Table S3). The TVE was significantly larger for negative compared with positive correlation at high intensity of warming and sediment deposition, whereas the opposite (not significant) pattern was observed at low treatment intensity (post-hoc contrasts, Table S3). Similarly to what observed for biomass, the experimental F I G U R E 3 Biomass simulations for different correlation scenarios of warming and sediment deposition. Panels show the response of biofilm biomass (as μg chl a cm −2 ) for different levels of warming, sediment deposition, and their spatial correlation in simulations. Black dots indicate the values of the two variables and their correlation chosen as treatment levels in the subsequent experimental test of simulation predictions. LMEMs indicated no differences between controls for artifacts and controls for the three response variables examined ( Figure S8, Table S5).
The correlation effect -i.e. the contribution of spatial correlation to nonlinear averaging -emerges when two stressors act non-additively, that is, when the effect of the concurrent change of two stressors is different from the sum of the effects of changing each stressor individually (Crain et al., 2008;Koussoroplis et al., 2017). In our simulations, the moderate correlation effect in the low-intensity condition originated from the mild antagonistic effect of warming and sediment deposition, where low levels F I G U R E 5 Testing predictions. Predicted vs. experimental values of biomass (μg chl a cm −2 ) as a function of the correlation and intensity of warming and sediment deposition. Yellow and dark-green filled circles indicate mean biofilm biomass under negative and positive correlation, respectively (n = 3). Gray circles refer to mean biomass of transects exposed to constant conditions (n = 3). Purple reversed-triangles indicate mean biomass of plots exposed to constant conditions derived from the first experiment (n = 3). Yellow and dark-green empty squares are the expected values of biofilm biomass obtained from simulations under low (∆T = +5°C and sediment deposition = 0.5 cm) and high (∆T = +15°C and sediment deposition = 1.5 cm) intensity, respectively. Error bars are non-parametric bootstrapped 95% confidence intervals. In particular, error bars of predictions incorporate the effect of error propagation stemming from uncertainty associated with the RS. of sediment deposition (0.5 cm) may have buffered the adverse effects of warming through nutrient release and reduction of desiccation stress (Figure 4a) (Larson & Sundbäck, 2012;McKew et al., 2011;Nakamoto et al., 2000). The correlation effect evanished with increasing sediment deposition, weakening the interaction between sediment and warming. Increasing sediment loads may have generated hypoxic conditions at the sediment-biofilm interface, such that the antagonistic buffering effect caused by the 0.5 cm layer of sediment switched into a negative effect operating additively with warming under the thicker 1.5 cm layer.
The correlation effect reversed under high warming and sediment deposition, with a positive spatial correlation of the two variables resulting in adverse effects on both the TVE and biofilm biomass compared with negative spatial correlation (Figures 4c and 5a,b).
Such inversion of the correlation effect could also reflect a switch from an antagonistic to a synergistic effect of warming and sediment deposition on biofilm biomass. However, caution must be paid in interpreting this inversion, as no previous experiments have directly characterized the nature (antagonistic or synergistic) and the magnitude effect of warming and sediment deposition on rocky intertidal biofilm.
The observed differences among correlation scenarios might also reflect a shift in species composition. Although physiological responses to stress may occur within a few days in microorganisms, longer periods (months) may be necessary for these effects to translate into compositional shifts (Schimel et al., 2007). In a parallel experiment, we found that changes in species composition in response to warming required at least 4 months to be detected (L. Rindi, unpublished data). Since the experiment presented here lasted 3 months, we are more inclined to believe that outcomes were driven more by physiological responses than shifts in species composition within the biofilm.
Although the direct manipulation of warming and sediment reproduced the effects of intensity and spatial correlation observed in the simulations, RS predictions may be affected by error propagation ( Figure 5). This likely reflected the uncertainty associated with estimating the RS from field data. However, results of the bootstrapping analysis showed that despite a moderate fit (R 2 = 0.29), our RS generated realistic predictions of response of biofilm to changes in intensity and spatial correlation of warming and sediment deposition.
Our study examined the effects of nonlinearities and correlation between variables at small spatial scales. Whether our RS could predict the TVE at larger spatial scales remains an open question.
Scale-transition theory uses Jensen's inequality (nonlinear averaging) and measures of environmental variances and correlations to extrapolate local ecological patterns (e.g., population dynamics within patches of habitat) to broader scales (e.g., regional population dynamics) Chesson et al., 2005;Melbourne & Chesson, 2006). The effect of Jensen's inequality increases with the degree of nonlinearities and with the amount of variance encountered when embracing larger spatial or temporal scales. Thermal performance curves are a typical example of nonlinear response functions that have been widely used to model population responses in fluctuating environments (Denny, 2017;Kingsolver & Woods, 2016;Koussoroplis & Wacker, 2016). Studies have shown how the shape of performance curves may change depending on the duration and history of exposure to stressful conditions, acclimatation, and ontogeny that may all affect the accuracy of temporal predictions from Jensen's inequality (Kingsolver & Woods, 2016;Kremer et al., 2018;Sinclair et al., 2016). Similar effects may occur when extrapolating across spatial scales. Due to the patchy nature of the rocky intertidal environment, biofilm assemblages might become increasingly different in terms of disturbance legacies, acclimatation, and overall response to warming and sediment deposition with increasing spatial scales, compromising the ability of the RS to predict patterns at larger scales. Although these caveats remain to be clarified, our study shows how consideration of nonlinear response functions and spatial correlation can help elucidating the influence of multiple processes at small spatial scales in a heterogeneous environment.
Focusing on local patterns and processes is important because small-scale variability is ubiquitous in nature, and most of the interactions between organisms and the surrounding environment occur at the scale of the microhabitat (Korell et al., 2021;Pincebourde et al., 2016;Potter et al., 2013). Yet, microhabitat heterogeneity can promote adaptation by buffering populations against adverse environmental conditions and small-scale processes can drive largescale patterns of community stability (Grman et al., 2010;Riddell et al., 2021). For example, ridges and depressions may generate various patterns of correlation among solar radiation, freezing, and moisture in tundra systems, providing favorable microclimates for the persistence, growth, and adaptation of dominant shrub species that are responsible for coarse-scale vegetation shifts (greening) in Arctic and Alpine ecosystems (Dobbert et al., 2021). Similar smallscale patterns of correlation between leading environmental variables are expected to occur in other terrestrial and aquatic systems where landscape features and topographic complexity promote finescale environmental heterogeneity (Deák et al., 2021). Assessing the generality of nonlinear and correlation effects as mechanisms shaping small-scale ecological patterns is important to better understand the compound effects of multiple processes and the ecological role of microclimates in changing environments.
Climate change projections for the 21st century involve modifications of the spatiotemporal patterns of climate variables (Gunderson et al., 2016;Hayashida et al., 2020;Young & Ribal, 2019). The spatial context in which organisms are embedded, such as landscape and microtopographic features, strongly filter and modify the climatechange signal (Pincebourde et al., 2016). The thermal mosaics originating from the topographic complexity of rocky intertidal shores are an example of a filtered signal (Helmuth et al., 2006). In the same vein, filtered signals may originate in a multifactorial context when topographic features generate small-scale patterns of correlation of environmental variables. As we have shown here, local patterns of correlation may play a crucial role in modulating the effect of climate extremes, ultimately influencing average biofilm biomass at the scale of the shore. Our analysis signals the need for researchers, resource managers, and policymakers aiming at predicting the impact of multiple stressors to account for current and possibly future spatiotemporal patterns of correlation among stressors. Black swans can occur in space as well as in time (Anderson & Ward, 2019), but the spatial context is often overlooked in the analysis of ecological extremes. Spatial correlation should be explicitly incorporated into multiple-stressors studies (Gunderson et al., 2016), risk-assessment framework (Côté et al., 2016;Goussen et al., 2016) and coupled environmental-physiological models (Pincebourde et al., 2016;Rezende et al., 2014). Our study shows how the inclusion of correlation between drivers can improve predictions from Jensen's inequality in real-world conditions. Further experimental work is needed to evaluate the generality of these findings in other ecosystems and over a wider range of environmental variables.

ACK N OWLED G M ENTS
We thank Chiara Ravaglioli for constructive comments to the manuscript and Caterina Mintrone for field and technical assistance.
The research was supported by the University of Pisa under projects PRA_2017_19 and PRA_2020_76.

CO N FLI C T O F I NTE R E S T
Authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are available on Figshare at the following https://doi. org/10.6084/m9.figshare.14447871.